E/ B decomposition of finite pixelized CMB maps 



Emory F. Bunn 
Physics Department, University of Richmond, Richmond, VA 23173 

Matias Zaldarriaga 

Physics Department, New York University, 4 Washington Place, New York, NY 10003 
Institute for Advanced Study, Einstein Drive, Princeton, New Jersey, NJ 08540 

Max Tegmark 

Dept. of Physics, Univ. of Pennsylvania, Philadelphia, PA 19104 

Angelica de Oliveira-Costa 
Dept. of Physics, Univ. of Pennsylvania, Philadelphia, PA 19104 



Separation of tlie E and B components of a microwave background polarization map or a weak 
lensing map is an essential step in extracting science from it, but when the map covers only part 
of the sky and/or is pixelized, this decomposition cannot be done perfectly. We present a method 
for decomposing an arbitrary sky map into a sum of three orthogonal components that we term 
"pure E" , "pure B" and "ambiguous". The fluctuations in the pure E and B maps are due only 
to the the E and B power spectra, respectively, whereas the source of those in the ambiguous map 
is completely indeterminate. This method is useful both for providing intuition for experimental 
design and for analyzing data sets in practice. We show how to find orthonormal bases for all three 
components in terms of bilaplacian eigenfunctions, thus providing a type of polarized signal-to-noise 
eigenmodes that simultaneously separate both angular scale and polarization type. The number of 
pure and ambiguous modes probing a characteristic angular scale 9 scales as the map area over 9^ 
and as the map boundary length over 9, respectively. This implies that fairly round maps (with 
short perimeters for a given area) will yield the most efficient i?/B-decomposition and also that 
the fraction of the information lost to ambiguous modes grows towards larger angular scales. For 
real-world data analysis, we present a simple matrix eigenvalue method for calculating nearly pure 
E and B modes in pixelized maps. We find that the dominant source of leakage between E and B 
is aliasing of small-scale power caused by the pixelization, essentially since derivatives are involved. 
This problem can be eliminated by heavily oversampling the map, but is exacerbated by the fact 
that the E power spectrum is expected to be much larger than the B power spectrum and by the 
extremely blue power spectrum that CMB polarization is expected to have. We found that a factor 
of 2 to 3 more pixels are needed in a polarization map to achieve the same level of contamination 
by aliased power than in a temperature map. Oversampling is therefore much more important for 
the polarized case than for the unpolarized case, which should be reflected in experimental design. 



I. INTRODUCTION 



Detecting polarization of the cosmic microwave background (CMB) radiation has become one of the main goals 
of the CMB community. Numerous experimental groups are currently searching for CMB polarization 0-|5j. CMB 
polarization can potentially offer a vast amount of information about our Universe. In general, polarization in 
very sensitive to the ionization history of the Universe. For example, on large scales it can provide insight into 
the way the Universe reionized 1^. On degree scales, once the temperature anisotropics are well measured, the 
predicted polarization can serve as a test of how and when recombination happened and could potentially lead to an 
important confirmation of the Big Bang model 0,^. Moreover, because the bulk of the polarization is produced at 
the last-scattering surface, it should exhibit no correlation on scales larger than about one degree unless there were 
super-horizon perturbations at decoupling. Polarization can thus become a good test of inflation Q . 

Most of the recent interest in polarization is based on its ability to provide evidence for a stochastic background of 
gravity waves. It has been shown that the polarization field on the sky can be decomposed into two parts, a scalar part 
usually called E and a pseudoscalar part usually called B [|lO|,fll . The pseudoscalar part cannot be created by density 



perturbations to linear order in perturbation theory. A detection of the B component on large scales would thus 



indicate the presence of a background of gravity waves, a prediction of inflationary models [y_2 13 1. Such a detection 
would determine the energy scale of inflation and could provide a stringent test of inflationary models p^ ] . On smaller 
scales, the B modes will most probably be dominated by secondary contributions produced after last scattering, the 
leading one being gravitational lensing p^ . A detection of these contributions could provide information about the 
distribution of matter all the way up to the last-scattering surface. There are many proposals for how to detect and 
use this effect p^-p^. In standard models, however, the B component is likely to be quite difficult to detect [p^-pl|. 

It is clear that a separation of the observed polarization into E and B parts is crucial to much of the CMB 
polarization scientific program. It has been realized, however, that real-world complications such as the finite size of 
the observed patch can significantly reduce our ability to do a clean separation between the two components: when 
using a quadratic estimator method for measuring the E and B power spectra, substantial "leakage" between the two 
was found on large angular scales . In it was shown that naive estimates of the sensitivity needed to detect the 
B component that ignore the such leakage can significantly underestimate the required sensitivity for an experiment 
aimed at detecting the B modes. In it was shown that in a finite patch, modes that are only E or only B can be 
constructed but that there are also ambiguous modes, modes that receive contributions to their power from both E 
and B. The construction of the modes was done for a round patch working in harmonic space. It was shown for each 
value of m there are two ambiguous modes. 

The issue of separating E and B has also generated interest in the field of weak gravitational lensing [ p6| - 
where the basic cosmological signal is expected to produce only an i?-pattern in cosmic shear maps, and the B-mode 
therefore serves as an important test for other signals due to intrinsic galaxy alignment or systematic errors. Although 
we do not discuss weak lensing explicitly in this paper, our results are relevant to that case as well since the lensing 
E/B problem is mathematically analogous. 

In this paper we revisit the issue of E and B mode separation, with two goals: to provide intuition for experimental 
design and for efficiently analyzing data sets in practice. We present a general derivation of the pure E, pure B and 
ambiguous modes in real space, and relate them to the eigenfunctions of the bilaplacian on a finite patch. We then 
introduce a way to obtain modes that are very nearly "pure" in a pixelized map by solving a generalized eigenvalue 
problem and discuss how this can be used to analyze real-world data sets. 

The paper is organized as follows. Section II establishes some notation and reviews the mathematics underlying 
the E/B decomposition of a polarization field. In Section III, we show how to decompose the space of all polarization 
fields on a finite patch of sky into pure E modes, pure B modes, and modes that are ambiguous with respect to the 
E/B decomposition. Section IV presents examples of this decomposition. In Section V, we present a method for 
finding (nearly) pure E and B modes numerically for pixelized maps by solving a generalized eigenvalue problem. 
Section VI presents examples. In Section VII we show that aliasing of small-scale power is the dominant source of 
"leakage" between the i^-modes and the i3-modes. We summarize our conclusions in Section VIII. 



II. E AND B MODES: NOTATION AND PRELIMINARIES 



In this section we will review the definition of E and B modes to introduce all the relevant notation. We will also 
give alt erna tive definitions of these modes which will help clarify how this decomposition works on finite patches of 
sky. In |II C] we discuss the small-angle approximation. Further details on properties of spin-two fields on the sphere 
and the E/B decomposition may be found in, e.g., ||23| , ^ and references therein. 



A. Spin two notation 



This section is rather technical. Since all intuitive aspects of our results can be understood in terms of t he m uch 
simpler formulae that apply in the flat-sky approximation, some readers may wish to skip straight to Section |ll C| and 
revisit this section as needed. 

The (linear) polarization of the CMB is described in terms of the Stokes parameters Q and U. The definition of Q 
and U depends on the coordinate system chosen. In this subsection we review definitions that are valid for the full 
sky, so we will use spherical coordinates to define Q and U. 

We will follow the notation of ||T^. The Stokes parameters can be combined to form a spin 2 {Q + iU) and a spin 
—2 (Q — iU) combination. In the full sky these combinations can be decomposed using spin-2 harmonics. 
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It is natural to introduce a scalar {E) and a pseudoscalar {B) field to describe polarization. The expansion 
coefficients of these two fields in (ordinary spin-0) spherical harmonics are 



+ a_2,im)/2 



i{0'2,lm — a-2,im)/2. 



(2) 



On the sphere, these two functions completely characterize the polarization field They are important physically 
because cosmological density perturbations cannot create B type polarization while gravitational waves can pO| , pl| . 
On small scales B polarization can be generated by lensing , and furthermore B may turn out to be a good monitor 
of foreground contamination, although at the moment nothing is known about how different foregrounds contribute 
to E or B. In terms of aE.im and asdm the Stokes parameters can be written as p5[| 



Q = - y^(Q£,iTi 
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-2^;m)/-^- These functions satisfy X^i^ = —Xij^m and 



'X2,i-m which together with a*^ = aE,i-m and a 



o.B,i-m make Q and U real quantities. 



The spin-2 harmonics in ec^uation (|1|) can be related to the usual spin-0 spherical harmonics by means of two first- 
order differential operators, the spin-raising {d ) and spin-lowering ) operators [0, which are defined in spherical 
coordinates by 
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where s is the spin of the function to which the operator is being applied. When applied to the spin- weighted spherical 
harmonics, these operators yield the following identities: 



d sYun = -[{1 + s)il - s + ^_ii^,„. 
In particular, the spin-0 and spin-2 harmonics are related as follows: 

2Yi^^[{l-2)\/{l + 2)\]^/^ddYi^, 

-2Y,„, = [ii-2y./ii + 2y.]'/'33Yi^. 

Another useful consequence of these relations is 

{l + 2)\ 



ddddYi,n = d333 Y, 



Yi^^{l + 2)il + l)l{l-l)Yi^, 
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or equivalently that when acting on spin-zero variables, 

33dd ^3333 =v2(v2 + 2), 

since corresponds to + 1) in spherical-harmonic space. 
Equations (|l]), (|) and (|) can be combined to obtain: 

Q + iU = 3 3 i^/jE + i^b) ; Q ~ tU = 3 3 i^/jE - ii'B) 

r^E = - - 2)1/(1 + 2)l]'^\E,lmYlm ; = - - 2)1/(1 + 2)\]'^^aB,lmYl„^. 
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Thus Q and U can be written in terms of second derivatives of the scalar and pseudoscalar "potentials" ipE and ipB, 
which are directly related to E and B. Equation ([l0| ) is analogous to the fact that a vector field can be written 
as a sum of a gradient and a curl component. The difference for spin-2 fields is that one can write them as second 
derivatives of the scalar and pseudoscalar potentials. 

We pause to note that the reason why E and B are the focus of attention instead of tpE and ipB is partly a matter 
of convention. Perhaps more importantly, E and B have the same power spectrum on small scales as the Stokes 



parameters, while the derivatives in equation (|lO|) imply that the power spectra of the Stokes parameters and those 
of i^E and differ by a factor [l - 2)!/(/ + 2)! - r*. 

To clarify the relation between all these quantities, we can think of weak gravitational lensing {e.g., [^,^). The 
shear variables are the analogues of the Stokes parameters, E is the analogue of the projected mass density, and 
is the analogue of the projected gravitational potential. 

We can use equations (|lO| ) and (||) to show that, 

V'(V' + 2)i;E ^[dHQ + ^U)+dd{Q~ iU)]/2 

V2(V2 + 2)Vb = i[dd {Q + iU)~dd{Q~ iU)]/2. (11) 

These equations show that we can take linear combinations of second derivatives of the Stokes parameters and obtain 
variables that depend only on E or on B. (In the flat-sky approximation, the left-hand sides of these equations are 
simply V'^E and V^B respectively. On the sphere, the relation is not so simple, but it is still true that the left-hand 
sides depend only on E and B respectively.) We will use this to project out the E and B contributions. 



B. Vector notation 



We can summarize the above results using a slightly different notation that will help clarify the analogy with vector 



D. = U ) (12) 



D^^if^^f -ff ^ ). (13) 



fields. We will use boldface to denote the polarization field written in the form of a vector, P = f/ ) ■ then 
define two second-order differential operators and D^;, 

2 \ -i{dd -dd) 

1 / i{d d ~ d 6 

2 V 66+66 

Equation ( |lo|) now becomes 

P = D£;fe+DBl/;B, (14) 

the analogue of the gradient/curl decomposition. Moreover, D^; and satisfy two important properties, 

D^j • Ds = D^j . Db = 0, (15) 
D^-D£ = Dtj.DB=v2(v2 + 2). (16) 

Equation (|l|) is the spin-two analogue of the familiar fact that V x V = 0. Substituting equation (|l0| ) into equa- 
tion (15) implies that if a polarization field on the sky has only i? as a source, it should satisfy • P = and if it 

is only due to a i? component is should satisfy • P = 0. 
In this vector notation, equation (^) can be written as 

P = — ^ O.EJm ^EJm + 0,B,lm B ,lnn 
Im 

Xl^lm \ V ( ^^2.1 

■E.l-n 



C. Small-angle approximation 

In this subsection, we present some formulas valid in the small-angle (flat-sky) approximation. When working in 
this limit, it is more natural to measure the Stokes parameters with respect to a Cartesian coordinate system (cc, y) 
instead of the usual polar coordinate axis. In the flat-sky approximation, the differential operators reduce to simply 
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Using the above expressions it is trivial to demonstrate that -Ds = -Ds = and that -Ds — -Ds = V''. 
In the flat-sky approximation, jV^j S> 1 (that is, only modes with eigenvalues much greater than one contribute 
significantly), so the V^CV^ + 2) operator in equation (^6|) has reduced to the bilaplacian V^. 

De and Ds are the spin-2 analogues of the familiar gradient and curl operators. Applyin g D g or Ds to a scalar 
field gives E and B fields that have vanishing "curl" and "gradient", respectively. Equations (20) and (21) show that 
Db = RDb, where the 2x2 matrix 
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(22) 



simply performs a rotation taking Q ^ —U and U Q. When drawing polarization fields as two-headed arrows with 
length {Q^ + U^Y^^ and angle tan^^(J7/Q)/2, this corresponds to rotating the polarization direction by 45° at each 
point . In other words, rotating the polarization directions of an i?-field by 45° gives a i?-field. 
The analogue of equation (n% is now given in terms of Fourier modes. 
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(23) 



In other words, the E/B decomposition becomes local in Fourier space: the polarization direction of the i?-component 
is parallel or perpendicular to k whereas that of the S-component makes a 45° angle with k. 



III. A NATURAL BASIS FOR POLARIZATION FIELDS 



On a manifold without boundary, any polarization field can be uniquely separated into an E part and a B part. 
But if there is a boundary (i.e., if only some subset of the sky has been observed), this decomposition is not unique. 
Let us first introduce some notation to clarify the problem. 

Polarization fields living on form a normed vector space with the inner product 

(P,P')= / P-P'dn, (24) 

and we say that two fields P and P' are orthogonal if (P, P') = 0. We refer to a polarization field P as 

• i? if it has vanishing cur 1, I.e., • P = 0, 

• S if it has vanishing divergence, i.e., • P = 0, 

• pure E if it is orthogonal to all i?-fields, and 

• pure B if it is orthogonal to all i?-fields. 

As long as fl is simply connected, which we shall assume throughout this paper, an equivalent definition of an E 
polarization field is one that can be derived from a potential ipE via P = 'De'>Pe- (And, of course, an analogous 
statement holds for B fields. As always, the analogy with the more familiar case of vector fields holds: any curl-free 
field is the gradient of a potential.) 

On the complete sky, every polarization field can be uniquely represented as a linear combination of an E field and 
a B field, and all E fields are perpendicular to all B fields. In other words, the space of all polarization fields is the 
direct sum of two orthogonal subspaces: the space of all E fields and the space of all B fields. [One way to prove 
these assertions is simply to use the E and B spherical harmonics defined in ( p7| ) as a basis.] In this case, there is no 
distinction between an E field and a "pure E'' field. 



But if only some subset of the sky has been observed, so that is a manifold with boundary, then this decomposition 
is not unique. One way to see this is to note that there are modes that satisfy both the iJ-mode and B-mode conditions 
simultaneously. When we split a polarization field into an E part and a B part, these "ambiguous" modes can go into 
either component. In order to make the E/B decomposition unique, we must first project out the ambiguous modes. 

In other words, the subspaces of all E modes and all B modes are no longer orthogonal: in fact, they overlap. To 
recapture orthogonality, we must restrict our attention to the pure E and B subspaces. To be specific, the space of 
pure E modes is the orthogonal complement of the space of all B modes, which includes both pure B modes and 
ambiguous modes. Similarly, the space of pure B modes is orthogonal to both the pure E modes and the ambiguous 
modes. In summary, we can represent the space of all polarization fields on f2 as a direct sum of three subspaces: 
pure E, pure B, and ambiguous. 

In this section, we show explicitly how to construct orthonormal bases of pure E modes, pure B modes, and 
ambiguous modes, so that we can unambiguously decompose any polarization field into these three components. In 
|p2| this construction was presented for a cap working in harmonic space. We here present the general formalism 
in real space. For simplicity, we work in the flat-sky approximation, although the construction works without this 
assumption. 

We first construct the ambiguous modes. An ambiguous mode P must be an i?-mode, so P = D^/ for some scalar 
field /. And it must also satisfy the B-mode condition: • P = 0. Combining these, we get 



So we can make a pair of ambiguous modes D^/ and DbI out of any function / that satisfies V^/ = 0. All such 
biharmonic functions are determined by their values and first derivatives on the boundary of the region, so it is 
straightforward to form a basis of them simply by choosing a basis for the set of scalar functions on the boundary. 

In the quest of separating the E and B contributions the ambiguous modes are not very useful, since we cannot 
know whether they are due to a cosmological E ot B signal. If we are willing to assume (on either observational or 
theoretical grounds) that E dominates over B on the angular scale of interest, then it may be sensible to assume that 
power found in the ambiguous modes is E power. This does enhance the accuracy with which the E power spectrum 
can be detected in a given data set [ pH . 

Of much more use are the "pure" E and B modes. We now give an explicit construction of these pure modes. 

Let the scalar field TpE generate a pure £^-mode GEi'E, and let Db'>Pb be any i?-mode (not necessarily pure). The 
requirement for a pure £'-mode is that these be orthogonal: 



If we use the explicit forms ( po|) and ( pi| ) for the differential operators and integrate by parts twice to move D e over 
to the DsV-B term, this reduces to a line integral around the boundary of fl. (After integrating by parts, the surface 
integral vanishes because it contains • D^V's, which is zero.) The line integral contains terms proportional to ipE 
and n • VipE- The conditions for a pure E mode are therefore 

1. -^B = on the boundary Oft. 

2. n • WipE = on the boundary dft. 

In other words, '0B must satisfy both Dirichlet and Neumann boundary conditions simultaneously. Fortunately, 
the bilaplacian operator has a complete set of eigenfunctions that satisfy these boundary conditions. To form an 
orthogonal basis of pure E modes, all we have to do is find a complete set of such eigenfunctions and apply the 
operator D^; to them. Similarly, if we apply D^, we will have an orthogonal set of pure _B-modes. The pure B modes 
can also be found by taking the pure E modes and rotating the polarization at each point by 45°. 

The boundary conditions turns out to have a simple geometrical interpretation: for a pure i?-mode the polarization 
on the boundary must be parallel or perpendicular to the boundary; for a pure i?-mode it must make a 45° angle 
with the boundary. 

The proof that these basis functions are orthogonal is similar to the more familiar situation with eigenfunctions 
of the Laplacian. Let / and g be eigenfunctions of with eigenvalues A and n, and let them satisfy the boundary 
conditions. Then D^;/ and D^jgr are two of our "pure i?" basis functions. Their inner product is 



= • Db/ = VV- 



(25) 




(26) 




(27) 



where we have integrated by parts twice and used the boundary condition on / to drop the boundary terms. Of course 
the same argument with / and g switched leads to the conclusion that the inner product is A times the integral of fg. 



If A 7^ /i, then the integral must therefore vanish, and if A = /i, we can take a hnear combination that orthogonahzes 
the two modes. We choose to normahze all modes P so that (P, P)=l. 

In conclusion, the pure i?-modes, pure i?-modes and ambiguous modes form a complete orthonormal basis for the 
space of all square-integrable {i.e., (P, P) < oo) polarization fields P in a sky region fi. We found that a polarization 
field is 

• pure E if it has vanishing curl and is parallel or perpendicular to the boundary, 

• pure B if it has vanishing divergence and makes a 45° angle with the boundary, and 

• ambiguous if it has vanishing divergence and curl. 

These conclusions apply not only to the eigenmodes that we have constructed but more generally, by linearity, to 
any field. This means that we can optionally decompose a polarization field P into its three components directly, 
without going through the step of expanding it in eigenmodes. The pure £^-componcnt Pb is obtained by solving the 
bi-Poisson equation = • P with Dirichlet and Neumann boundary conditions and computing P^; = 'De4'e- 

The pure _B-componcnt Ps is obtained analogously, and the ambiguous component P? is simply the remainder, i.e., 
P.=P-Pe-Pb 

IV. WORKED EXAMPLES I 

In this section we illustrate the above construction for two worked examples: a disk in the flat-sky approximation 
and a spherical cap. 

A. Disk 

Suppose that the observed region is a disk of radius R with R ^ 1 radian, so that the flat-sky approximation is 
appropriate. 

We begin with the ambiguous modes. We want to find functions / with V^f = 0. Assume a separable solution 
/(r, 0) = F(r)e™*. We know that V^(V^/) = 0, so V^/ must be a harmonic function. The most general solution is 
V^/ cx r™e"""^. Solving this equation for /, we get two independent solutions: 

fi-'^)^[ln.%^.^,- (28) 

As we expected, there are in general two solutions per m (since there are two conditions we wish to impose on the 
boundary). Each solution yields two ambiguous modes, D^/ and D^/, which turn out to be just rotations of each 
other. 

In the case to = 0, though, these two solutions do not yield any ambiguous modes, as TDeI = Ds/ = 0. The 
same is true for the first of the two to 1 solutions, so there is only one pair of ambiguous modes with to = 1. This 
counting of modes agrees with p^ . 

We now proceed to find the pure E- and pure _B-modes. One way to construct eigenfunctions of is to take 

f = ax + a-x, (29) 

where a\ is an eigenfunction of with eigenvalue A and a-\ has eigenvalue —A. (These two are obviously degenerate 
eigenvalues of V'' with eigenvalue A^ , so we can take linear combinations of them. Of course there are no well-behaved 
eigenfunctions of with positive eigenvalue over an entire manifold, but there are over a manifold with boundary.) 




FIG. 1. Ambiguous modes of a disk. From top to bottom, m = 1, 2, 3. 



Once again we apply separation of variables in polar coordinates. The angular dependence is e""'^. Then for any 
positive k, the Bessel function Jm{kr) has eigenvalue —k^ and the modified Bessel function Im{kr) has eigenvalue fc^, 
so we can take our eigenfunctions of to be 

fmk{r, <t>) = (aJmikr) + blm{kr))e'""^ . (30) 

The boundary conditions tell us that 

a _ I^kR) _ I^jkR) 

b JmikR) JLikR)' ^ ' 

So there will be solutions for all values of k that satisfy J^/Jm = I'ml^m- These roots can be computed numerically. 
For large n, a good approximation for the nth root with azimuthal quantum number m is 

(Tfl \ 
n + - j . (32) 

The figures show the first few modes of each type. As noted above, there arc no ambiguous modes with m = 0, 
one pair of ambiguous modes with |m| = 1, and two for each |m| > 1. Only one of each pair is shown; the other is 
found by rotating the whole pattern. Similarly, for each pure E and pure B mode, a linearly independent mode can 
be obtained by rotating the page. 





FIG. 2. The first two m = pure E (left) and pure B (right) modes for a disk. 



FIG. 4. Same as Fig. | with m = 2. 
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If our data covered the entire plane, we would construct a basis out of only the ordinary Bessel functions J, 
excluding the modified Bessel functions /„j. In the limit kR — > c», therefore, we expect the contribution from /,„ to 
be small, and indeed this is the case. The function grows exponentially for large argument, so in order to satisfy 
the Dirichlet boundary condition the coefhcient b in equation (^) must be small. For a mode with kR ^ 1, therefore, 
the ordinary Bessel function dominates except near the boundary. In this limit, the modified Bessel function takes 
over in a small region near the boundary to "flatten out" the mode and make it satisfy the Neumann boundary 
condition. 

It is worth noting that all modes except those with rn = require that both Q and U be measured in the patch. 
Modes with m — Q depend only on Q for the pure E modes and only on U for the pure B modes (with Q and U 
defined with respect to the polar coordinates). 



B. Spherical cap 

This construction can be adapted to give the basis functions for a spherical cap without recourse to the flat-sky 
approximation. In this case, the functions we are looking for are eigenfunctions of the operator 



V2(V2 + 2) = (V^ + if - 1. 



(33) 



The ambiguous modes will therefore be eigenfunctions of the Laplacian with eigenvalues and —2. These eigen- 
functions can be written in terms of associated Legendre functions Pim as 



•^^'"•^(^''^)-\Pi,„(cos0)e™* (34) 

for any integer m. 

The associated Legendre function Pim is well-behaved over the entire sphere as long as \m\ < I, so there appear to 
be four singularity-free solutions over the entire sphere. These are mapped to zero by D^; and T>b, though, so they 
do not give ambiguous modes. This is of course as it should be: there are no ambiguous modes over the entire sphere. 

If, however, the region of interest is a spherical cap < Q, then we permit functions that have singularities 
outside ri. In that case, there is one nontrivial ambiguous mode with m = ±1, namely Po±ie^"^, and two for every 
m with |m| > 1. The I = modes can be written explicitly as 

„ / „^ / sin 6 , , 

Pom(cos0)= — . (35) 



1 + COS 6* 



The / = 1 modes are not so simple. The first one is 



(cos0-t-2)sin^g 

Pi2 cos e) = — — —- , 36 

(1 + cos 0)'' 

and the remainder can be computed from recurrence relations. 

We can also construct the pure E and B modes from the associated Legendre functions. Suppose we fix the 
azimuthal quantum number m and look for eigenfunctions of V^(V^ + 2) with eigenvalue k. We can construct one 
by taking a linear combination of two associated Legendre functions P\^m and P\_rm where A± are the two roots of 

(-A±(A±+1) + 1)2 = A: + 1. (37) 

The left-hand side is the eigenvalue of (V^ -f 1)^; recall that the eigenvalue associated with Pi„i is + 1), and 
compare this equation to equation (^3[). 

Just as in the case of the disk, there will be a discrete set of A;'s for which a linear combination of these two functions 
can satisfy both boundary conditions. 



V. PIXELIZED MAPS 



In this section we study the decomposition of polarization in finite pixelized maps. One possibility would be to 
search for eigenfunctions of a discretized version of the bilaplacian operator. On scales much larger than the pixel 
scale, we would expect to recover modes that are approximately the same as those found above. The orthogonality 
of pure E and B modes would not be expected to be perfect in the discretized case, but on reasonably large scales it 
should be close. The main drawback of this approach is that by construction it explicitly assumes that both Q and 
U are measured at each pixel, so we would like to generalize the approach preserving its spirit and power. 

We will adopt a different method in which a complete set of E, B, and ambiguous modes can all be found at once 
by solving a single eigenvalue problem. With this approach, we can find a basis of modes that approximate the pure 
E and B modes very well (except for modes with frequencies close to the Nyquist frequency, where problems may be 
expected to arise no matter what approach one adopts). 

We assume that we have a map of a finite portion of the sky composed of N pixels. In each pixel we could have 
measured both Q and U; however, it is possible that in some or all of them only one combination of the Stokes 
parameters was measured. We will denote the vector of measured Stokes parameters P, which will have dimension 
less than or equal to 2N . In terms of the E and B modes of the full sky, the vector P is given by: 

P = — ^(flfi,/™ ^E,lm + asdm ^B,lm)- (38) 
Im 

Ideally we want to find the pixelized analogues of the pure E, pure B and ambiguous modes. A pure E mode, 
which we will denote e, should satisfy YB,im • e = for all Im. A pure B mode b satisfies YE^hn • b = for all Im. It 
is clear that one cannot find a solution to these sets of equations, i.e. to find such an e or b, since in general, we are 
trying to satisfy more equations than we have components of the P vector. In practice the number of constraints we 



need to satisfy is set by the angular resolution of the experiment, which determines the maximum / mode that has 
any appreciable power. Thus the difficulty of finding pure e or pure b modes will increase as the distance between 
pixels gets larger compared to the angular resolution of the experiment. Moreover, we also expect that the number of 
pure E and pure B modes will decrease as the fraction of pixels where only one of the Stokes parameters is measured 
increases. 



A pure E mode should satisfy 



or equivalently 



J2CBi\\YB,im-e\\^ = , (39) 



Im 



e*Be = ; B = J^CbiY B,i,nYl,i^, (40) 



Im 

for any choice of power spectrum Cbi- An analogous statement clearly holds for pure B modes: 

b*Eb = ; E = ^ Cb;Y£^,,^Y^ (41) 

Im 

for any Cei- The matrices E and B give the contribution to the power in each mode from the E and B components. 
In order to find candidate E and B modes numerically, we must choose a particular power spectrum; we will choose 
C(E.B)i/'^'^ — {I ~ 2)!/(Z + 2)! X W^, where is the window function that describes the beam smearing. We 
will motivate this choice in the next subsection: in practice we found it to work extremely well, making mixing 
between modes extremely small and almost perfectly recovering the modes we obtained in the previous section with 
the bilaplacian. 

Our aim is to construct a basis of vectors that span all the space but are ordered by their relative contributions 
from E and B modes. In principle, we would like to find the generalized eigenvectors of something like Ee = A^Be. 
A problem arises, however: we know that B has a null space (the space of pure E modes). So we regularize the 
problem by introducing a matrix N = ct^I, with I the identity matrix and a very small constant. We then solve 

(E + N)e = A£;(B + N)e. (42) 

If we choose small enough, the matrix E + N is essentially equal to to E in the subspace of pure E modes and is 
proportional to the identity matrix in the subspace of pure B modes. The converse holds for B + N. As a consequence, 
the eigenvectors with large A^; will be very close to pure E modes. Furthermore, with our choice of power spectra, 
the eigenvectors will automatically separate in scale with larger scale modes having a larger eigenvalue. 
There is an equivalent equation for B modes, 

(B + N)b = AB(E + N)b, (43) 

but any mode e satisfying equation ( |4^ ) also satisfies equation ( ^ ) with As = 1/Xe- 

We can derive simple and useful properties of the eigenvalues and eigenvectors if we assume that at every pixel in 
the map we have both Q and U. We consider the simple transformation where we rotate the polarization at every 
pixel by 45° {i.e. Q —U and U — > Q). We denote this transformation R45. It is represented by a block diagonal 
matrix 

(R45)., = S., ( 1 V ) ' (4^^) 
where i,j label pixels. The matrices E and B satisfy R45ER45 = B and R45BR45 = E. Moreover R45R45 = I. By 



substitution into equation (|4^), it is straightforward to prove that the vector e' — R45e also solves the eigenvalue 
equation but with eigenvalue 1/A^;. We conclude that if at every pixel we have measured both Q and U, modes that 
solve equation ( ^2|) come in pairs with eigenvalues Xe and 1/A^;. One member of the pair is preferentially E and the 
other preferentially B. 

In the next section we will present numerical examples to gain intuition on how the eigenvalue problem works. First 
we will motivate our choice of Ci spectra. 



A. Relation to bilaplacian formalism 



To find the relation between our eigenvalue and bilaplacian formalisms, we start by considering a vector satisfying 
the eigenvalue equation, 

(E + cr2l)e = AE(B + cr^I)e. (45) 
Since D^B = D^Eb = 0, multiplying equation ( ^5| ) by and yields two scalar equations, 

Dj5Ee = a2(A£-l)Dtje, 

D^Be = cr2(A£i - l)D|je. (46) 

We now proceed to show that with our choice of spectra Cei = (^ ~ 2)!/(Z + 2)! the modes constructed using our 
bilaplacian formalism solve equation (|6|). We take 

V'^iV'' + 2)i,E = \^E- (47) 

and assume that ■0b satisfies both Dirichlet and Neumann boundary conditions. We can use the completeness relation 
for spherical harmonics and our choice of spectra to write 

A-V2(v2 + 2)^^(0)- j de'Y,Yim{e)YC^{e')^E{e') 

Im 

= I d9'J2cEi^l±^Yi^ie)Y;^ie')^Ei9') (48) 

We can use the fact that D^Y^,;™ = [{I + 2)1/(1 - 2)l]'^/'^Yi,n to get 

f dO' CEi-DlYE.ira{d)Yl,^{e') BEi^E{e') - A" ^^(V^ + 2)V£(0) (49) 



Im 



where we have integrated by parts using the boundary conditions satisfied hy ipE- Finally we can factorize the 
bilaplacian operator • D^; = V^(V^ + 2) and use our definitions e — DEi^E and the E matrix to get 

D|;Ee = A^^D^e. (50) 

Thus if we identify A^^ ~ cr^i^E — 1), Db'0-E satisfies the first of equations (^). The second equation in (^) is 
trivially satisfied because being e a pure E mode it follows that both Be — and D^e — 0. 

We have just shown that modes constructed using the bilaplacian formalism solve equation ( ^6|) rather than (|4^). 
This means that the vector e = T^e^Je actually satisfies 

(E + a^I)e ^ X{B + aH)e + a, (51) 

where a has to be an ambiguous mode (because it has to give zero when acted upon by both and D^). The 
easiest way to understand what is happening is to look at the structure of the E and B matrices in the basis of the 
eigenfunctions of the bilaplacian. If we call a one of the basis vector in the ambiguous space and contract equation 
(|l|) with it we find, 

a*Ee = a • a, (52) 

where we have also used the fact that e was a pure E mode. Thus the reason why there is an extra ambiguous mode 
in eq. (51) is that the E matrix can have non-zero elements mixing the pure E and ambiguous subspaces. In other 
words our two formalisms are identical when restricted to the pure E and B subspaces but differ in the ambiguous 
subspace. 

In practice we will find that the modes calculated by solving the bilaplacian equation and the generalized eigenvalue 
problem are almost identical. This can be understood by looking at equation ( ^ ) and realizing that in most cases we 
will be able to achieve very good separation, i.e., A ^ 1. This implies that one only needs to add a very tiny amount 
of ambiguous modes to e in equation (51) to "correct it" and make a zero (because A is so large). This is especially 
so because under most circumstances the matrix elements of both E and B in the subspace of ambiguous modes are 
comparable. 



VI. WORKED EXAMPLES II 



We begin by revisiting the cap example we solved in the continuous case. We start by assuming that every pixel 
has both Q and U . We consider a fiducial experiment with a 0.2 degrees FWHM for the beam angular resolution. 
The patch observed has a radius of 3.8 degrees and contains 351 pixels (the spacing between pixels in both the radial 
and the tangential directions was set to 0.2 degrees as well). 




FIG. 5. E/B eigenvalues for a cap. In the top panel we show A_e as a function of mode number. On the bottom we show 
the eigenvalues as a function of both and l^^ as defined in equation (^5|). We took = A x 10~^, a factor 10~^ smaller 
than the zero lag correlation function. 

Figure ^ shows the eigenvalues we obtained. As expected, the eigenvectors come in pairs with eigenvalues A^; and 
1/Xe- The eigenvectors with very small eigenvalues correspond to pure B modes and those with very large ones to 
pure E modes. The particular values of the eigenvalues should not be given much importance as they depend on value 
of the regularizing constant a^. What is important is that the large eigenvalues show the good degree of separation 
that we have achieved. 

There is also a concentration of modes at A = 1. These modes have two origins. First, modes on small scales, 
where our small cr^I regularization dominates over the E and B matrices, will have A^; = 1. Second, as we discussed 
in the previous sections, there are larger-scale ambiguous modes that receive contributions from both E and B. Our 
method is unable to separate between both types because they have the same eigenvalues. 



FIG. 6. Examples of modes for a cap. We show the nine modes with lowest A_b. 



In Figure we show the first nine eigenvectors, corres pondi ng to the lowest nine eigenvalues. One immediately 
recognizes in this set the pure B modes discussed is section IV A . The first eigenvector corresponds to the lowest-order 
TO = mode. The next two are the lowest to = 1 modes, which differ only by a rotation. Then come the second 
TO = mode, then the lowest to = 2 modes, then the lowest to = 3 modes, and finally the second to = 1 mode. The 
best nine E modes, corresponding to the largest nine eigenvalues, are simply equal to the ones plotted in figure ^ but 
with each polarization "vector" rotated by 45°. 



FIG. 7. On the left (right) we show the best two B type eigenvectors for a 32 x 32 (16 x 64) pixel patch. 



Our method for finding modes can be used for any shape of sky patch. In figure [7| we show the first two modes of 
a square patch 32 x 32 pixels on a side. Comparing with figure ^ it is clear that they are essentially the same modes 
as the two first modes for the cap. We also show the first two modes in a patch 16 x 64 pixels on a side. We have 
also checked that these modes for the rectangle can be derived from the bilapacian formalism. 

To understand where the ordering of modes in figure ^ is coming from, i.e., why the modes appear in that order in 
the figure, we will introduce window functions for each mode. We define 



where we have introduced pi = l{l + l)C//27r. Using the window functions we can define an effective / for each mode, 
the average / calculated using the window function as a weight. Specifically, we can define quantities 




(53) 



(54) 
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that give the average I for the E and B contribution to a given mode. 



1 ^ 

0.8 
.0.6 
0.4 
0.2 



1 

0.8 
.0.6 t- 
0.4 r 
0.2 j- 

P- 

1 

0.8 
.0.6 
0.4 

0.2 

1 

0.8 
.0.6 t- 
0.4 r 
0.2 r 
r 




H 1 — I I I I I 1 1 1 1 — I I I I I I 



''I I ^ 



m=1 




I'l 



- 



H 1 — I I I I I I I 1 1 — I I I I I I I 1 1 — I I I I I 



-+- 



m=2 



I I I I I I 1 H 




II 
I I 
/ I 
I I 
' 1/ 



hH \ — [ 



m=3 



I - 



...VIA, 



-I I I I ■ ■ ■ ■ ' I I I I ■ ■ ■ ■ 



J I I l_l_L 



1 



10 



100 



1000 



FIG. 8. Window functions for the modes shown in figure |^. Sohd fines show Wi^ and dasfied lines sfiow Wf' , normalized to 
unit peak height. 



Figure |^ shows the window functions for the eigenvectors that were plotted in figure Jq. Note that the window 
functions are well localized in I and each has a clear peak. Moreover, the modes in figure ffare ordered in increasing 
order of Z^. Both of these are a consequence of our choice of power spectra. 



VII. ALIASING 

The windows Wf' in figure ^ can be used to determine where the leakage between E and B is coming from. The 
dotted line in all the panels gives an estimate of the Nyquist frequency in the map. The conclusion is clear: the 
contaminating power is aliased power. For power that is aliased one cannot distinguish E from B. The remedy for 
this is to increase the sampling in the map so as to further suppress the aliased power. 

The bottom panel of figure |5| shows the eigenvalues we obtained for the cap but as a function of Z^ and l^g. We see 
that the modes with large values of A^; have a low and a large Z^, indicating that most of the contamination is 
coming from aliasing. The opposite is true for modes with low A^;. We also see that some of the modes with A^; = 1 
receive contributions from large scales, an indication that these are truly ambiguous modes. 

We can understand our results intuitively by considering a simple toy model, closely following the treatment in 
We work in the small-angle limit and compute the Fourier components of the observed polarization field assuming 
they were observed over a square patch of size L. Using equation (^3|) we obtain 



Eici) 
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where we have defined the window function 



^ ' \L ) smkyAe/2' ^ ' 

with A9 the separation between pixels. The Nyquist wavenumber is /cNyq — tt/A^?. 

The naive way to recover the E and B components would be to combine the Fourier coefRcients P as one would 
do if the patch were infinite. For B for example, we would compute 

B{k) = - sin 2</> Q(k) + cos 2<?!) U{k) (58) 

and then estimate the power spectrum by taking the square of these variables. 
In terms of the real E and B, our _B-estimate can be written as 

S(k) = / -^iy(k-q)[sin2a £;(q)+cos2a B(q)], (59) 

J V'T^r 

where cos a = k • (\/kq. This estimate has contributions from both E and B. Only when 14^ is a delta- function such 
that a = do we avoid mixing. The E contributions arise because of two effects: the finite size of the sky patch and 
the pixelization (causing aliasing). The effect of the finite patch size manifests itself as a finite width of the peaks of 
the window function, while the effect of aliasing is that W has several peaks. 

We have shown how to construct modes that avoid contamination due to the finite patch size. These modes are 
not Fourier modes. In what follows we want to show that even in the limit L — s- cx), there is still mixing due to 
pixelization. If we take this limit the window function becomes a sum of delta functions centered at q = 2(m, n)fcNyq, 
where m and n are integers. 

To consider a concrete example, we calculate the ratio of power in E and _B in a mode with wave vector ko produced 
by an initial field that only had E modes with a power spectrum Ci and assuming an infinite but pixelized sky map. 
All the wave vectors k.y = kg + 2{i, j)k^y^ will contribute to this mode. We get 

(^(kp)) ^ E.,-sin^2(0,,--0o)Q.^./C;„ 

(^(ko)> Ey-cos2 2(0,,-</)o)Q.,/Q,' ^ ^ 

where Uj = |kij| and Iq = |ko|. Equation ( |60| ) shows that all the aliased modes contribute to B contamination because 
in general these modes do not have — (pio) = mi:. 

It is important to note that the aliased power is suppressed by the beam. As i and j become larger, the magnitude 
of the power on those scales decreases because Ci is proportional to Wi , the beam window function. For example, 
if we consider a mode with wavevector along the positive x axis, the aliased mode with the smallest possible beam 
suppression has a power suppressed by a factor exp[— 2fcNyq(fcNyq — IkoD/c^], where ui, is the Gaussian width of the 
beam (ab is related to the full width half max of the beam (FWHM) by Ob = ^^8 ln(2)/FWHM). For fixed ko, the 
suppression can be made as large as one wants by increasing fcNyq, that is by increasing the sampling of the map. If we 
want the beam to produce a suppression factor S we need to choose 2k^^^{l — |ko|/fcNyq)/(7^] = ln(5'), or equivalently 

FWHM/A6' « 0.5^1n(5')/(l- |ko|/fcNyq). 

A point worth noting about aliasing is that the power spectrum of the polarization is a rapidly growing function of I 
and that the E power spectrum is expected to be much larger than the B one. Figure ^ shows the power spectrum for 
E and B type polarization from gravity waves in a ACDM model. The temperature spectra were COBE normalized 
and the tensor component was assumed to be 10% of the temperature anisotropics on COBE scales. The sharp 
increase in power between E and B partially compensates the smearing by the beam. To give a rough feeling of what 
sampling is needed to avoid aliasing we could assume that we want the aliased power to be a factor of 100 smaller than 
the power we want to measure. For a temperature map that would correspond to a suppression factor S ^ 10^ while 
for polarization we would need S ~ 10^ which means that the ratio FWIIM/A0 has to be a factor \/2.5 larger for 
polarization than for temperature, or equivalently that one needs a factor of 2.5 more pixels to obtain the same level 
of contamination. We conclude that one has to be particularly careful about aliasing when dealing with polarization 
maps if one wants to obtain a clean separation between E and B. 




FIG. 9. Polarization power spectra in a ACDM model. The solid curve is for E produced by density perturbations and the 
dashed curve for the B component produced by tensor modes. The anisotropics were COBE normalized and it was assumed 
that the tensor component was 10% of the anisotropies on these scales. 



The effect of aliasing can be decreased by increasing the sampling of the map. It should be noted, however, that the 
presence of holes, bad pixels or pixels with only one measured Stokes parameter in the map will have a similar effect. 
We illustrate this by considering a toy example. We artificially increase the noise variance (the diagonal elements of 
N) for a fraction of the pixels chosen at random. Figure |lO| shows the window functions for the first nine modes in an 
example where 20% of the measured Stokes parameter were assumed to have the large noise. For comparison, we also 
show the original window functions. On large scales, the modes look essentially the same as the ones plotted in figure 
^. The effect of the missing pixels is very noticeable in the E window function, the one that quantifies the leakage. As 
might have been expected, the level of contamination coming from modes of frequency around the Nyquist frequency 
is greatly increased. 




FIG. 10. Effect of measuring only one Stokes parameter on 20% of the pixels chosen at random. The panels show the window 
functions for the first modes with m = 0,1,2,3 when both Stokes parameters are measured (solid lines) and when 20% are 
missing (dashed lines). 



VIII. DISCUSSION 



We have developed a formalism for measuring the E and B components of polarized CMB maps or weak lensing 
maps given the real-world complications of finite sky coverage and pixelization. 

We have shown that by expanding a map in a particular basis, obtained by differentiating bilaplacian eigenfunctions, 
it can be decomposed as a sum of three orthogonal components that we term pure E, pure B and ambiguous. The 
pure £^-component is orthogonal to all B-modes and are therefore guaranteed to be caused by an E signal (on the 
uncut sky), and conversely for the pure B-component. The ambiguous component is the derivative of a biharmonic 
function, and the original map contains no information about whether it is due to E- or _B-signal in the uncut sky. 
We also derived a discrete analogue of these results, applicable to pixelized sky maps. Our results are useful both for 
providing intuition for survey design and for analyzing data sets in practice. 



A. Implications for survey design 



To maximize our ability to separate E and B, we clearly want to minimize the fraction of modes that are ambiguous. 
We found that the ambiguous modes are specified along the boundary of the map rather than in the two-dimensional 
interior. This means that the number of pure and ambiguous modes probing a characteristic angular scale 6 scales as 
the map area over 6^ and as the map boundary length over 9, respectively. It is therefore best to minimize the ratio 
of circumference to area, i.e., to make the patch as round as possible. 



Almost all pure modes (all except the ones with m = for the spherical cap example) are a combination of both Q 
and U Stokes parameters, so to achieve unambiguous E/B separation, one needs to measure both, with comparable 
sensitivity throughout the map. 

With pixelized maps, we found that aliasing of small-scale power was a serious problem. Although it can in principle 
be eliminated by heavily oversampling the map, the required oversampling is greater than for the unpolarized case, 
both because derivatives are involved and because CMB polarization is expected to have an extremely blue power 
spectrum. This has important implications for, e.g., the Planck satellite, where bandwidth constraints on the telemetry 
have been mentioned as reasons to reduce the oversampling. It is crucial to bear in mind that the usual Nyquist 
rule-of-thumb that applies to unpolarized maps may be insufficient for realizing the full scientific potential of Planck's 
CMB polarization measurements because one needs roughly a factor of 2 to 3 more pixels in a polarization map to 
achieve the same level of contamination by aliased power. 



B. Implications for data analysis 

In [ pO| , it was shown how a quadratic estimator method could produce uncorrelated measurements of the E and 
B power spectra from real-world data sets with arbitrary sky coverage, pixelization and noise properties, and this 
method has been apphed to both the POLAR and PIQUE [| data. The one annoying problem with this method 
was that it gave E/B leakage. Our present results allow us to understand and eliminate this problem. 

We now know that leakage is caused by the ambiguous modes. The above-mentioned scaling tells us that the 
fraction of modes probing a given angular scale I that are ambiguous scales as l^-^, in good agreement with the 

asymptotic behavior empirically found in |20| . Although presented a technique for removing most of the leakage, 
we now know how to remove it completely: by eliminating the ambiguous modes. 

In practice, the way to do this is to compute two projection matrices He and Ub that project onto the subspaces 
given by the eigenvectors e of equation with A^; > A* and A_b < 1/A,, respectively, for some large eigenvalue cutoff 
A*, say A* — 100. The three maps n^P, n^P and [I — IIe — HbIP will then be approximately the pure E, pure B 
and ambiguous components of the original map P, which can be directly used for visual inspection, cross-correlation 
with other maps and systematic error tests. To measure the E and B power spectra, one compresses the original data 
vector P into two shorter ones P^; and Pb by expanding it into the above-mentioned pure E and pure B eigenvectors, 
respectively. Since this is a mere matrix multiplication, the corresponding noise and signal covariance matrices (which 
the quadratic estimation method takes as input) are trivially computed as well. These two data vectors will each have 
less than half the length of P. Since the time required by the quadratic estimator method scales as n^, the final E and 
B power spectrum calculations are therefore about an order of magnitude faster than in the original |^o|] approach. 

It should be noted that the ambiguous modes are not useless in all circumstances. If it has been established that E 
dominates over B (as is expected theoretically) by observing the pure modes, then it is safe to assume that most of the 
power in the ambiguous modes is E power as well. In this case, the ambiguous modes can be used to reduce the errors 
on estimates of the E power spectrum. This could be particularly useful when attempting to constrain reionization 
with .E-power on the very largest angular scales attainable with a galaxy-cut all-sky map, where a substantial fraction 
of the modes will be ambiguous. 
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